Learning Goals

Lab Description

We will work with COVID data downloaded from the New York Times. The dataset consists of COVID-19 cases and deaths in each US state during the course of the COVID epidemic.

The objective of this lab is to explore relationships between cases, deaths, and population sizes of US states, and plot data to demonstrate this

Steps

I. Reading and processing the New York Times (NYT) state-level COVID-19 data

0. Install and load libraries

1. Read in the data

# Combine data by state
cv_states_readin <- as.data.frame(data.table::fread("https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv") )


state_pops <- as.data.frame(data.table::fread("https://raw.githubusercontent.com/COVID19Tracking/associated-data/master/us_census_data/us_census_2018_population_estimates_states.csv"))

state_pops$abb <- state_pops$state
state_pops$state <- state_pops$state_name
state_pops$state_name <- NULL

cv_states <- merge(cv_states_readin, state_pops, by="state")

2. Look at the data

  • Inspect the dimensions, head, and tail of the data
  • Inspect the structure of each variables. Are they in the correct format?
dim(cv_states)
## [1] 39114     9
head(cv_states)
##     state       date fips   cases deaths geo_id population pop_density abb
## 1 Alabama 2021-08-21    1  659750  12000      1    4887871    96.50939  AL
## 2 Alabama 2021-09-17    1  764839  13048      1    4887871    96.50939  AL
## 3 Alabama 2020-07-24    1   76005   1438      1    4887871    96.50939  AL
## 4 Alabama 2021-09-03    1  714860  12394      1    4887871    96.50939  AL
## 5 Alabama 2020-03-18    1      51      0      1    4887871    96.50939  AL
## 6 Alabama 2022-02-09    1 1254032  17452      1    4887871    96.50939  AL
tail(cv_states)
##         state       date fips cases deaths geo_id population pop_density abb
## 39109 Wyoming 2021-04-09   56 56873    701     56     577737    5.950611  WY
## 39110 Wyoming 2021-09-18   56 83958    918     56     577737    5.950611  WY
## 39111 Wyoming 2020-06-01   56   910     17     56     577737    5.950611  WY
## 39112 Wyoming 2021-04-20   56 57456    705     56     577737    5.950611  WY
## 39113 Wyoming 2020-04-11   56   343      0     56     577737    5.950611  WY
## 39114 Wyoming 2021-01-06   56 45890    464     56     577737    5.950611  WY
str(cv_states)
## 'data.frame':    39114 obs. of  9 variables:
##  $ state      : chr  "Alabama" "Alabama" "Alabama" "Alabama" ...
##  $ date       : IDate, format: "2021-08-21" "2021-09-17" ...
##  $ fips       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ cases      : int  659750 764839 76005 714860 51 1254032 159169 932250 3953 104786 ...
##  $ deaths     : int  12000 13048 1438 12394 0 17452 2558 16503 114 1882 ...
##  $ geo_id     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ population : int  4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 ...
##  $ pop_density: num  96.5 96.5 96.5 96.5 96.5 ...
##  $ abb        : chr  "AL" "AL" "AL" "AL" ...

3. Format the data

  • Make date into a date variable
  • Make state into a factor variable
  • Order the data first by state, second by date
  • Confirm the variables are now correctly formatted
  • Inspect the range values for each variable. What is the date range? The range of cases and deaths?
cv_states$date <- as.Date(cv_states$date, format="%Y-%m-%d")


state_list <- unique(cv_states$state)
cv_states$state <- factor(cv_states$state, levels = state_list)
abb_list <- unique(cv_states$abb)
cv_states$abb <- factor(cv_states$abb, levels = abb_list)


cv_states = cv_states[order(cv_states$state, cv_states$date),]


str(cv_states)
## 'data.frame':    39114 obs. of  9 variables:
##  $ state      : Factor w/ 52 levels "Alabama","Alaska",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ date       : Date, format: "2020-03-13" "2020-03-14" ...
##  $ fips       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ cases      : int  6 12 23 29 39 51 78 106 131 157 ...
##  $ deaths     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ geo_id     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ population : int  4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 ...
##  $ pop_density: num  96.5 96.5 96.5 96.5 96.5 ...
##  $ abb        : Factor w/ 52 levels "AL","AK","AZ",..: 1 1 1 1 1 1 1 1 1 1 ...
head(cv_states)
##       state       date fips cases deaths geo_id population pop_density abb
## 571 Alabama 2020-03-13    1     6      0      1    4887871    96.50939  AL
## 380 Alabama 2020-03-14    1    12      0      1    4887871    96.50939  AL
## 431 Alabama 2020-03-15    1    23      0      1    4887871    96.50939  AL
## 33  Alabama 2020-03-16    1    29      0      1    4887871    96.50939  AL
## 168 Alabama 2020-03-17    1    39      0      1    4887871    96.50939  AL
## 5   Alabama 2020-03-18    1    51      0      1    4887871    96.50939  AL
tail(cv_states)
##         state       date fips  cases deaths geo_id population pop_density abb
## 38547 Wyoming 2022-03-18   56 155907   1769     56     577737    5.950611  WY
## 39053 Wyoming 2022-03-19   56 155907   1769     56     577737    5.950611  WY
## 38906 Wyoming 2022-03-20   56 155907   1769     56     577737    5.950611  WY
## 38959 Wyoming 2022-03-21   56 155907   1769     56     577737    5.950611  WY
## 39020 Wyoming 2022-03-22   56 155988   1783     56     577737    5.950611  WY
## 38587 Wyoming 2022-03-23   56 155988   1783     56     577737    5.950611  WY
head(cv_states)
##       state       date fips cases deaths geo_id population pop_density abb
## 571 Alabama 2020-03-13    1     6      0      1    4887871    96.50939  AL
## 380 Alabama 2020-03-14    1    12      0      1    4887871    96.50939  AL
## 431 Alabama 2020-03-15    1    23      0      1    4887871    96.50939  AL
## 33  Alabama 2020-03-16    1    29      0      1    4887871    96.50939  AL
## 168 Alabama 2020-03-17    1    39      0      1    4887871    96.50939  AL
## 5   Alabama 2020-03-18    1    51      0      1    4887871    96.50939  AL
summary(cv_states)
##            state            date                 fips           cases        
##  Washington   :  793   Min.   :2020-01-21   Min.   : 1.00   Min.   :      1  
##  Illinois     :  790   1st Qu.:2020-09-05   1st Qu.:16.00   1st Qu.:  49810  
##  California   :  789   Median :2021-03-12   Median :29.00   Median : 209010  
##  Arizona      :  788   Mean   :2021-03-12   Mean   :29.78   Mean   : 536429  
##  Massachusetts:  782   3rd Qu.:2021-09-16   3rd Qu.:44.00   3rd Qu.: 665336  
##  Wisconsin    :  778   Max.   :2022-03-23   Max.   :72.00   Max.   :9079164  
##  (Other)      :34394                                                         
##      deaths            geo_id        population        pop_density       
##  Min.   :    0.0   Min.   : 1.00   Min.   :  577737   Min.   :    1.292  
##  1st Qu.:  884.2   1st Qu.:16.00   1st Qu.: 1805832   1st Qu.:   43.659  
##  Median : 3525.5   Median :29.00   Median : 4468402   Median :  107.860  
##  Mean   : 8812.6   Mean   :29.78   Mean   : 6419723   Mean   :  422.717  
##  3rd Qu.:10688.8   3rd Qu.:44.00   3rd Qu.: 7535591   3rd Qu.:  229.511  
##  Max.   :88461.0   Max.   :72.00   Max.   :39557045   Max.   :11490.120  
##                                                       NA's   :741        
##       abb       
##  WA     :  793  
##  IL     :  790  
##  CA     :  789  
##  AZ     :  788  
##  MA     :  782  
##  WI     :  778  
##  (Other):34394
min(cv_states$date)
## [1] "2020-01-21"
max(cv_states$date)
## [1] "2022-03-23"

4. Add new_cases and new_deaths and correct outliers

  • Add variables for new cases, new_cases, and new deaths, new_deaths:
    • Hint: You can set new_cases equal to the difference between cases on date i and date i-1, starting on date i=2
  • Filter to dates after July 1, 2021
for (i in 1:length(state_list)) {
  cv_subset = subset(cv_states, state == state_list[i])
  cv_subset = cv_subset[order(cv_subset$date),]

  # add starting level for new cases and deaths
  cv_subset$new_cases = cv_subset$cases[1]
  cv_subset$new_deaths = cv_subset$deaths[1]


  for (j in 2:nrow(cv_subset)) {
    cv_subset$new_cases[j] = cv_subset$cases[j] - cv_subset$cases[j-1] 
    cv_subset$new_deaths[j] = cv_subset$deaths[j] - cv_subset$deaths[j-1]
  }
  
  cv_states$new_cases[cv_states$state==state_list[i]] = cv_subset$new_cases
  cv_states$new_deaths[cv_states$state==state_list[i]] = cv_subset$new_deaths
}

# Focus on recent dates
cv_states <- cv_states %>% dplyr::filter(date >= "2021-07-01")
  • Use ggplotly for EDA: See if there are outliers or values that don’t make sense for new_cases and new_deaths. Which states and which dates have strange values?
p1 <- ggplot(cv_states,aes(x=date,y=new_cases, colour=state)) +
  geom_line() +  
  ylab("New Cases") +
  xlab("Date") +
  labs(title="New COVID Cases Overtime by State") +
  theme_minimal()

ggplotly(p1)
p2 <- ggplot(cv_states, aes(x=date,y=new_deaths,colour=state)) +
  geom_line() +
  ylab("New Deaths") +
  xlab("Date") +
  labs(title="New COVID Deaths Overtime by State") +
  theme_minimal()

ggplotly(p2)
  • Correct outliers: Set negative values for new_cases or new_deaths to 0

It was strange that there are negative deaths and cases, which were results of redefining terms by state authorities.

  • Inspect data again interactively
cv_states$new_cases[cv_states$new_cases < 0] = 0
cv_states$new_deaths[cv_states$new_deaths < 0] = 0

p3 <- ggplot(cv_states, aes(x=date, y=new_deaths, color = state)) +
  geom_line() +
  geom_point(size=0.5,alpha=0.5) +
  ylab("New Deaths") +
  xlab("Date") +
  labs(title="New COVID Deaths Overtime by State (outliers removed)")+
  theme_minimal()

ggplotly(p3)

5. Add additional variables

  • Add population-normalized (by 100,000) variables for each variable type (rounded to 1 decimal place). Make sure the variables you calculate are in the correct format (numeric). You can use the following variable names:

    • per100k = cases per 100,000 population
    • newper100k= new cases per 100,000
    • deathsper100k = deaths per 100,000
    • newdeathsper100k = new deaths per 100,000
  • Add a “naive CFR” variable representing deaths / cases on each date for each state

  • Create a dataframe representing values on the most recent date, cv_states_today

cv_states$per100k =  as.numeric(format(round(cv_states$cases/(cv_states$population/100000),1),nsmall=1))
cv_states$newper100k =  as.numeric(format(round(cv_states$new_cases/(cv_states$population/100000),1),nsmall=1))
cv_states$deathsper100k =  as.numeric(format(round(cv_states$deaths/(cv_states$population/100000),1),nsmall=1))
cv_states$newdeathsper100k =  as.numeric(format(round(cv_states$new_deaths/(cv_states$population/100000),1),nsmall=1))

cv_states = cv_states %>% mutate(naive_CFR = round((deaths*100/cases),2))

max_date <- max(cv_states$date)
cv_states_today = cv_states %>% filter(date==as.Date(max_date))

II. Scatterplots

6. Explore scatterplots using plot_ly()

  • Create a scatterplot using plot_ly() representing pop_density vs. various variables (e.g. cases, per100k, deaths, deathsper100k) for each state on most recent date (cv_states_today)
    • Color points by state and size points by state population
    • Use hover to identify any outliers.
# Restrict sizes to a range, specify how we want the markers to be
cv_states_today %>%
  plot_ly(x=~pop_density, y=~cases,type="scatter", mode="markers", color=~state,
          size=~population, sizes=c(5,70),marker=list(sizemode='diameter'),
          opacity=0.5) %>%
    layout(title = 'Cases by Population Density and State',  xaxis = list(title = 'Population Density'), 
         yaxis = list(title = 'Cases'), legend = list(title=list(text='<b> State </b>')))
  • Remove those outliers and replot.
# Remove the outlier not because it's an error, but just to make observing the other states' data easier
cv_states_today %>%
  filter(state != "District of Columbia") %>%
  plot_ly(x=~pop_density, y=~cases, color=~state, type="scatter", mode="markers", 
          size=~population, sizes=c(5,70), marker=list(sizemode="diameter", opacity=0.5)) %>%
  layout(title = 'Cases by Population Density and State',  xaxis = list(title = 'Population Density'), 
         yaxis = list(title = 'Cases'), legend = list(title=list(text='<b> State </b>')))
  • Choose one plot. For this plot:
  • Add hoverinfo specifying the state name, cases per 100k, and deaths per 100k, similarly to how we did this in the lecture notes
  • Add layout information to title the chart and the axes
  • Enable hovermode = "compare"
# Paste0 puts spaces in between the categories
cv_states_today %>%
  filter(state != "District of Columbia") %>%
  plot_ly(x=~pop_density, y=~cases, color=~state, type="scatter", mode="markers", 
          size=~population, sizes=c(5,70), 
          marker=list(sizemode="diameter", opacity=0.5), hover_info="text",
                      text=~paste(
                        paste0("State: ", state), 
                        paste0("Cases per 100k: ", per100k),
                        paste("Death per 100k:", deathsper100k),
                      sep="<br>")) %>%
  layout(title = 'Population-normalized Cases per 100k',  
         xaxis = list(title = 'Population Density'), 
         yaxis = list(title = 'Cases per 100k'), 
         hovermode="compare",
         legend = list(title=list(text='<b> State </b>')))

7. Explore scatterplot trend interactively using ggplotly()

  • For pop_density vs. newdeathsper100k create a chart with the same variables using ggplotly()
  • Explore the pattern between \(x\) and \(y\)
  • Explain what you see. Do you think pop_density correlates with newdeathsper100k?
p4 <- cv_states_today %>%
  filter(state != "District of Columbia") %>%
  ggplot(aes(x=pop_density, y=newdeathsper100k, colour = state, size=population)) +
  geom_point() +
  ylab("New Deaths per 100k") +
  xlab("Population Density") +
  labs(title="Population Density versus New Deaths per 100k") +
  theme_minimal()

ggplotly(p4)

8. Multiple line chart

  • Create a line chart of the naive_CFR for all states over time using plot_ly()
    • Use the zoom and pan tools to inspect the naive_CFR for the states that had an increase in September. How have they changed over time?
  • Create one more line chart, for Florida only, which shows new_cases and new_deaths together in one plot. Hint: use add_layer()
    • Use hoverinfo to “eyeball” the approximate peak of deaths and peak of cases. What is the time delay between the peak of cases and the peak of deaths?
cv_states %>%
plot_ly(x=~date, y=~naive_CFR, color=~state, type="scatter", mode="lines") %>%
  layout(title = 'Naive CFR Overtime',  
         xaxis = list(title = 'Date'), 
         yaxis = list(title = 'Naive CFR'), 
         legend = list(title=list(text='<b> State </b>')))
cv_states %>%
  filter(state=="Florida") %>%
  plot_ly(x=~date, y=~new_cases, type="scatter", mode = "lines", 
          name = "Cases") %>%
  add_lines(x=~date, y=~new_deaths, type="scatter", mode="lines",
            name="Deaths") %>%
    layout(title = 'Florida\'s COVID Cases and Deaths Overtime',  
         xaxis = list(title = 'Date'), 
         yaxis = list(title = 'New Cases'), 
         legend = list(title=list(text='<b> Type </b>')))

9. Heatmaps

Create a heatmap to visualize new_cases for each state on each date greater than July 1st, 2021 - Start by mapping selected features in the dataframe into a matrix using the tidyr package function pivot_wider(), naming the rows and columns, as done in the lecture notes - Use plot_ly() to create a heatmap out of this matrix. Which states stand out?

cv_states_mat <- cv_states %>%
  select(state, date, new_cases) %>%
  filter(date > "2021-07-01")

# Change pivot so that we have date as a column
cv_states_mat2 <- as.data.frame(pivot_wider(cv_states_mat, names_from=state, 
                                           values_from=new_cases))
# Here, colnames is not a variable so we don't need to write a function
cv_states_mat2 <- cv_states_mat2 %>%
  column_to_rownames("date") %>%
  as.matrix()

# Z value for color intensity
plot_ly(x=colnames(cv_states_mat2), y=rownames(cv_states_mat2), z=~cv_states_mat2, type="heatmap") %>%
      layout(title = 'Heatmap for New Cases by State Overtime',  
         xaxis = list(title = 'State'), 
         yaxis = list(title = 'Date'))

Deliverables

Lab 10b questions 1-2, lab 11 questions 0-10. Upload html or pdf for both lab Rmd’s to quercus.